{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5c2f94e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd \n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import bisect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a1ac26d",
   "metadata": {},
   "source": [
    "# Model\n",
    "\n",
    "The model is\n",
    "$$\n",
    "\\max_{m_i} \\frac{((1-\\tau)y_i-\\mu pm_i)^{1-\\sigma}}{1-\\sigma} + \\phi (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "The FOC is \n",
    "$$\n",
    "\\mu p ((1-\\tau)y_i- \\mu pm_i)^{-\\sigma} = \\phi\\alpha_0 \\exp(-\\alpha_1 - \\alpha_0 m_i)\n",
    "$$\n",
    "This gives a solution $m_i = {\\cal M}(y_i,p|\\sigma, \\phi, \\alpha_0, \\alpha_1, \\mu)$\n",
    "\n",
    "Health insurance system\n",
    "$$\n",
    "\\tau \\sum_i \\omega_i y_i = (1-\\mu)p \\sum_i \\omega_i m_i~~~\\Rightarrow~~~\\tau = (1-\\mu)s\n",
    "$$\n",
    "where $s=pm/y$ with $m=\\sum_i \\omega_i m_i$ and $y=\\sum_i \\omega_i y_i$\n",
    "\n",
    "\n",
    "### US\n",
    "For the US, we have $p=1$ and $\\sum_i \\omega_i y_i \\equiv y = 1$. Therefore, the moments restrictions  \n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i} = \\sum_i \\omega_i m_i\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i \\omega_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "$$\n",
    "\\Psi_3 = \\frac{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\max}))}{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\min}))}\n",
    "$$\n",
    "allows to identify $\\{\\phi,\\alpha_0,\\alpha_1\\}$ with $\\tau = (1-\\mu) \\Psi_1$\n",
    "\n",
    "\n",
    "### Europe\n",
    "\n",
    "For Europe, we keep $\\{\\phi,\\alpha_0,\\alpha_1\\}$ found in the previous step, and the moments restrictions\n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i \\omega_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "allows to identify $\\{p,\\alpha_1\\}$, given that a calibration of $\\sum_i \\omega_i y_i \\equiv y$ and $\\tau = (1-\\mu) \\Psi_1$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e15e6936",
   "metadata": {},
   "source": [
    "We start by programming functions for finding optimal $m$ from the first-order condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "01fd442f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def foc(m,sigma_c,phi,alphas,y,p,mu,mom):\n",
    "\treturn mu*p*( (1-(1-mu)*mom)*y - mu*p*m )**(-sigma_c) - alphas[0]*phi*np.exp(-alphas[1]-alphas[0]*m)\n",
    "\n",
    "def u(m,sigma_c,phi,alphas,y,p,mu,mom):\n",
    "\treturn ( ((1-(1-mu)*mom)*y - mu*p*m)**(1-sigma_c) )/(1-sigma_c) + phi*(1-np.exp(-alphas[1] - alphas[0]*m))\n",
    "\n",
    "def optm(sigma_c,phi,alphas,y,p,mu,mom):\n",
    "\ttry:\n",
    "\t\tm = bisect(foc,0.0,y*0.99,args=(sigma_c,phi,alphas,y,p,mu,mom))\n",
    "\texcept:\n",
    "\t\tup  = u(0.99*y,sigma_c,phi,alphas,y,p,mu,mom)\n",
    "\t\tlow = u(0.0,sigma_c,phi,alphas,y,p,mu,mom)\n",
    "\t\tif up > low:\n",
    "\t\t\tm = 0.99*y\n",
    "\t\telse :\n",
    "\t\t\tm = 0.0\n",
    "\treturn m\n",
    "\n",
    "def prod(alphas,m):\n",
    "\treturn 1 - np.exp(-alphas[1] - alphas[0]*m)\n",
    "\n",
    "def hopt(sigma_c,phi,alphas,y,p):\n",
    "\tm = optm(sigma_c,phi,alphas,y,p,mu,mom)\n",
    "\treturn prod(alphas,m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87cca3d1",
   "metadata": {},
   "source": [
    "These are the parameters that we end up getting, but we check here that the function works. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "72da1031",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_c   = 2\n",
    "p_us      = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1273ad26",
   "metadata": {},
   "source": [
    "We load up the parameters of the income process and use that to compute 4 median incomes within quartiles. The stationary variance of the income process of each country is "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8a9b58b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigs = [0.256367, 0.094643, 0.214538, 0.237439, 0.169834, 0.094643, 0.270357, 0.486429]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21276dcb",
   "metadata": {},
   "source": [
    "For Europe, we take a mean. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "54c54890",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.486429, 0.19111728571428574)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sig_us = sigs[-1]\n",
    "sig_eu = np.mean(sigs[0:7])\n",
    "sig_us,sig_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5096d3a0",
   "metadata": {},
   "source": [
    "We program up the income distribution, making sure it has a mean of 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "38b7dad4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.37922262, 0.67735473, 1.05644178, 1.88698087]), 1.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qs    = [0.125,0.375,0.625,0.875]\n",
    "yy_us = [norm(0,np.sqrt(sig_us)).ppf(q) for q in qs]\n",
    "yy_us = np.exp(yy_us)\n",
    "yy_us = yy_us/np.mean(yy_us)\n",
    "yy_us,np.mean(yy_us)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ce6034",
   "metadata": {},
   "source": [
    "We do the same for Europe but scale down everything to have an average income of 0.78"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5f4c9032",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.44109972, 0.63452033, 0.83837729, 1.20600267]), 0.78)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_eu  = 0.78\n",
    "qs    = [0.125,0.375,0.625,0.875]\n",
    "yy_eu = [norm(0,np.sqrt(sig_eu)).ppf(q) for q in qs]\n",
    "yy_eu = np.exp(yy_eu)\n",
    "yy_eu = yy_eu/np.mean(yy_eu)\n",
    "yy_eu = y_eu*yy_eu\n",
    "yy_eu,np.mean(yy_eu)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e279f29",
   "metadata": {},
   "source": [
    "We use the following moments for the U.S. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5bafd2b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_us    = 0.15\n",
    "h_us    = 0.89\n",
    "grad_us = 1.27\n",
    "mu_us   = 0.136"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "388ef5ca",
   "metadata": {},
   "source": [
    "We use the following moments for Europe (we do not use the gradient but still compare):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b139b8ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_eu    = 0.096\n",
    "h_eu    = 0.93\n",
    "grad_eu = 1.06\n",
    "mu_eu   = 0.155"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4ece99cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "no_social_security = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "14de71ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "if no_social_security:\n",
    "    mu_us = 1\n",
    "    mu_eu = 1\n",
    "else:\n",
    "    mu_us   = 0.136     \n",
    "    mu_eu   = 0.155"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e52353c8",
   "metadata": {},
   "source": [
    "# Calibration\n",
    "\n",
    "We start the exercise with finding U.S. parameters that fit the moments. We are solving a set of equations for three unknowns, $\\alpha_0$, $\\alpha_{1,US}$ and $\\phi$. We use the methods of moments to minimize the squared distance of moments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a9f55003",
   "metadata": {},
   "outputs": [],
   "source": [
    "def criterion_us(theta,sigma_c,alphas,y,p,mu,moms):\n",
    "\talphas[0] = theta[0]\n",
    "\talphas[1] = theta[1]\n",
    "\tphi = theta[2]\n",
    "\tms = [optm(sigma_c,phi,alphas,i,p,mu,moms[0]) for i in y]\n",
    "\tmoms_s = np.zeros(3)\n",
    "\tmoms_s[0] = p*np.sum(ms)/np.sum(y)\n",
    "\tmoms_s[1] = np.mean([prod(alphas,m) for m in ms])\n",
    "\tmoms_s[2] = prod(alphas,ms[3])/prod(alphas,ms[0])\n",
    "\treturn (moms[0] - moms_s[0])**2 + (moms[1] - moms_s[1])**2 + (moms[2] - moms_s[2])**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b05de664",
   "metadata": {},
   "source": [
    "These are the estimates reported in the text:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f419d91b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([7.19595858, 1.47802119]), 2.9541541221896157)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alphas_us    = np.zeros(2)\n",
    "opt_us       = minimize(criterion_us,[5.0,1.5,2],args=(sigma_c,[0.0,0.0],yy_us,p_us,mu_us,[s_us,h_us,grad_us]))\n",
    "alphas_us[0] = opt_us.x[0]\n",
    "alphas_us[1] = opt_us.x[1]\n",
    "phi          = opt_us.x[2]\n",
    "alphas_us,phi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2c30303",
   "metadata": {},
   "source": [
    "We can check that we match moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "02e9c4af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.15000552312921223, 0.8900057537200122, 1.270002390515131)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ms_us = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us]\n",
    "np.mean(ms_us), np.mean([prod(alphas_us,m) for m in ms_us]), prod(alphas_us,ms_us[3])/prod(alphas_us,ms_us[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "64e69769",
   "metadata": {},
   "outputs": [],
   "source": [
    "hs_us = [prod(alphas_us,m) for m in ms_us]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eacecdf",
   "metadata": {},
   "source": [
    "In Europe, we seek to estimate $p_{EU}$ and $\\alpha_{1,EU}$ using $(s_{EU},h_{EU})$ as moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bcf770d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def criterion_eu(theta,sigma_c,phi,alphas,y,mu,moms):\n",
    "\talphas[1] = theta[0]\n",
    "\tp = theta[1]\n",
    "\tms = [optm(sigma_c,phi,alphas,i,p,mu,moms[0]) for i in y]\n",
    "\tmoms_s = np.zeros(2)\n",
    "\tmoms_s[0] = p*np.mean(ms)/np.mean(y)\n",
    "\tmoms_s[1] = np.mean([prod(alphas,m) for m in ms])\n",
    "\treturn (moms[0] - moms_s[0])**2 + (moms[1] - moms_s[1])**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d05ed9b",
   "metadata": {},
   "source": [
    "These are the estimates reported in the paper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4936d7b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.867110640650491, 0.5407133717867146)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alphas_eu    = np.zeros(2)\n",
    "alphas_eu[0] = alphas_us[0]\n",
    "opt_eu = minimize(criterion_eu,[2.0,0.8],args=(sigma_c,phi,alphas_eu,yy_eu,mu_eu,[s_eu,h_eu]))\n",
    "alphas_eu[1] = opt_eu.x[0]\n",
    "p_eu = opt_eu.x[1]\n",
    "alphas_eu[1],p_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea21023f",
   "metadata": {},
   "source": [
    "We can check that we fit the targeted moments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "cb7573da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.09599366642858816, 0.9299734154871621)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ms_eu = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu]\n",
    "p_eu*np.mean(ms_eu)/np.mean(yy_eu), np.mean([prod(alphas_eu,m) for m in ms_eu])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b72224d5",
   "metadata": {},
   "source": [
    "The gradient is lower than in the U.S. but slighly larger than the target (1.06)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "7da1f53a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1318813651014679"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod(alphas_eu,ms_eu[3])/prod(alphas_eu,ms_eu[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1525ccb7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD4CAYAAADrRI2NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAl4klEQVR4nO3dd3iUVdrH8e9JQkISQggk9A4JTQhiKCJFBBFsWBADiiv2XSxr2RXd8rrL6mJbxbIq776uZaWJIFixUxUIkBAIkIQeaiiBQAgpc94/ZmQjIgyQ5Jny+1xXLjOZ5yT3PIYfh3vuecZYaxERkeAR4nQBIiJSvRT8IiJBRsEvIhJkFPwiIkFGwS8iEmTCnC7gZOLj423Lli2dLkNExG8sX758r7U2wZtjfTL4W7ZsSVpamtNliIj4DWPMFm+PVatHRCTIKPhFRIKMgl9EJMgo+EVEgoyCX0QkyCj4RUSCjIJfRCTIKPhFRBy27/AxZqdv57XvNlTLz/PJF3CJiASysnIX6dsKmJedz7zsfDK3H8RaaBRbkzv7tiIstGr35Ap+EZFqsKPgKPM9Qb8wdy+FxWWEhhjOb1aHhwYl0b9dAuc1jiUkxFR5LQp+EZEqUFxazrLN+5m33h32OXsOA+5d/RWdG9E/KYHebeOJjaxR7bUp+EVEKoG1lk17jxxv3/ywcR/FpS7Cw0Lo2aouN3ZvRv+kBNrWr4UxVb+rPxUFv4jIWTp8rIzFuXuZl53P/Jx8tu0/CkDr+GhSuzenf7sEerWqR2R4qMOV/pSCX0TES9ZasnYeYn72XuZl7yFt8wHKXJbo8FB6t43n7n5t6J+UQLO6UU6XekoKfhGRU9h/pIQFOfnMz97L/Jx88guPAdCxUW3u7NeafokJXNAijvAw/5mOV/CLiFRQVu4iI6+AednuFs6qvAKshbioGvRNTKBfUgL9EuOpX7um06WeNQW/iAS9XQeLj49aLsjJ51BxGSEGujarw28HukctOzeJJbQaRi2rg4JfRILOsbJy0jYfcE/grM9n/e5CABrWrsmQ8xrSP6k+fdrGExtV/aOW1UHBLyJBYXOFUcvvN+zjaGk54aEhdG8Vx/UXtKd/Un2SGjg/alkdFPwiEpCOHCvj+w37jo9abtlXBEDLelGMSGnqHrVsXY+o8OCLweB7xCLiW/ashfTJ0G4otOh91t/GWsu6XYXuoM/OZ9nm/ZSWW6LCQ+ndph539GlFv6QEWtSLrsTi/ZOCX0Sq35G9kDkDMibDzgwwoVCrwRkHf0FRCQty9h5/YnaPZ9SyfcMYbuvTiv6JCVzQMo6IMN96AZXTFPwiUj3KjkH255AxFXK+AFcZNOwCQybAecOhVsJpv0W5y5KRV3A86DO2FeCyEBtZg76J8fRLSqB/UgIN/HjUsjoo+EWk6lgL25e7WzmrP4DiAvfOvtevIXkkNOh02m+x+1DxT65qWVBUijGQ3LQO912SSP92CSQ3rRMwo5bVQcEvIpWvYBusmube3e/LgbCa0P5Kd9i3vhhCfzl6SspcpG3Zf3zUct0u96hlQkwEgzo0oH9SAn3axhMXHV5NDybwKPhFpHIcOwxr50DGFNi0ALDQvDdcdD90HAY1Y39x6dZ9RczL3sO87HwWb9hHUUk5NUINKS3qMm5oe/onJdC+YUxQjFpWBwW/iJw9VzlsXgDpU9yhX1oEcS3h4nHQ5Uao2+qky4pKyvhh477j16rf7Bm1bFY3kuu7NaV/UgIXtqlHdIQiqirorIrImcvPdu/sV02DQ9shojZ0vsHdymneC07YmVtryd59+PiuftmmA5SUu4isEcqFbeox5iL3qGXLelHa1VcDBb+IeKdov/sJ2owp7idsTQi0GQiDx0O7y6FG5E8OP1hUysJc9+WL52fvZdehYgDaNYjhV71b0D+pPikt46hZQ6OW1U3BLyK/rKwEcr90T+VkzwVXKdTvBIP/Bp1HQEyD44eWuyyZ2w8yb737lbIrtx7AZaF2zTD6JMbTP8l9ZctGsZGn+IFSHRT8IvJT1sKOle6JnNUzoGgfRCdAj7sgORUadTl+6J7CYvd16j1XtTzgGbXs0rQO9w5oe3zUMizUf65VHwwU/CLidmjHf0cw89dBaLi7hdN1FLS5BEJrUFLmYsXGfcdHLbN2HgIgvlYEA9rXp39SAn0TE6irUUufpuAXCWYlR2DdJ+5WzsbvAAvNesKVL0CnayEyjm37i5i3bId71DJ3L0dKygkLMVzQIo7fD2lH/6QEOjSsTYheQOU3FPwiwcblgi2L3Dv7rA+h5DDENod+v4PkVI7GtOSHTfuY9+UO5mdnsHHvEQCaxkVyzflNjo9axtQMzGvVBwOvgt8YMwSYCIQC/7LWTjjh/jjgTaANUAzcZq1d7blvM1AIlANl1tqUSqteRLy3b4N7IidjGhzcCuG1oOM12ORUciO7MC9nH/M+zGfJpmxKylxEhIVwYZt6jL6wBf2SEmgdH61RywBx2uA3xoQCrwKXAnnAMmPMHGttVoXDHgfSrbXXGmPae44fWOH+AdbavZVYt4h44+gBWDPL/QKrvKXuEczWF1PU73EWhPTk242HmT81nx0HFwKQWL8Wt/RyB32PVnU1ahmgvNnx9wByrbUbAYwxU4FhQMXg7wj8HcBau84Y09IY08Bau7uyCxaR0ygvhQ3fuPv26z+D8mPYhPbs6v4Yn5m+fLbVsOKDAspda4mJcI9a3jfQPWrZpI5GLYOBN8HfBNhW4XYe0POEYzKA64CFxpgeQAugKbAbsMAXxhgLvGGtnXSyH2KMuQu4C6B58+Zn8hhEBGDnKnffPnM6HMnHFVmPjc1v4EPbl8lb67J/WylQQJemsfy6fxv6t0uga7M61NCoZdDxJvhP1tSzJ9yeAEw0xqQDmcBKoMxz30XW2h3GmPrAl8aYddba+T/7hu6/ECYBpKSknPj9ReRkCne7gz5jKuxejSukBrl1LmJqSF/ezU+k9EAY9aLD6d/OfZ36PonxxNeKcLpqcZg3wZ8HNKtwuymwo+IB1tpDwBgA4372Z5PnA2vtDs9/9xhjZuFuHf0s+EXES6VHYf2nkD4Fu+FrjHWxMaID79nb+aCoB4XFtbmgeRy/7eYO+46NNGopP+VN8C8DEo0xrYDtQCowquIBxpg6QJG1tgS4A5hvrT1kjIkGQqy1hZ7PBwN/rcwHIBIUrIWtP1C2cjKsmUVYaSF7TD3eL72SmeV9Ka7Zln7JCUxISqB323rU1qilnMJpg99aW2aMuReYi3uc801r7RpjzD2e+18HOgDvGGPKcT/pe7tneQNglmcELAyYbK39vPIfhkhgsvs3sf/7d6mxehq1j+ZxzEbwuasHs+mPadmXfu0a8EZSPG0SamnUUrxmrPW9dnpKSopNS0tzugwRRxQW7GPzgsnUWvc+rY5k4LKGxa6OLIwehG1/NRd2aE7PVvWIDNeopfyXMWa5t6+T0it3RRzmclmytu9n09JPiMv5gJSji+hsStloGzOz7u2Y5BGkdOlCn7pRTpcqAULBL+KAfYePsSBnL9mrltBw8ywuc83nPFNAoalFVsOrCb/gZtp160/rMO3qpfIp+EWqQVm5i5XbCpifnU/6uhza7vmc60MWcE3IZsoJZXej/hzscTOxXa6kW5jGLaVqKfhFqsiOgqPMz3a/p+yS3J30LFnK9aELeCA0g7CwcoriO+NKeZrQzsNpHB3vdLkSRBT8IpVsUe5exn+cxbpdhzjf5DI6cjHPhy4mKrwQV62GhCTfC8kjiarfwelSJUgp+EUqSXFpOc98vp65i5ZyW+2lXF93IXWKtmBDIjEdroLkVEJaXwwh6tuLsxT8IpVg3eYdfDTtdQYd/pI/18yCEqBRH+j6e0yHq6FmbadLFDlOwS9ytlzluDbOZ8OXk2i+62t+Z45RFNsCuv8ButwIcS2crlDkpBT8Imcqfz1kTKEsfRphh3dQ30aRFjuYrlfeQ+3Ei0CvoBUfp+AX8UbJEVj5nvsdrHaswGVCWWST+dA1gouuGM31PdrokgniNxT8IqezeRHMHgsHNlFevzOzE8by1LbzaNGiJS+M6ErzenpFrfgXBb/ILyk5Al//FZa8AXEtWD3oP9y5IIr8wmP8dnAi9/RvQ5jexET8kIJf5GS2LIYPfwMHNlHW/S7+UZ7KPz/eReuEUGb+pjddmtZxukKRs6bgF6mopMizy38d6jRn61XTuWtBJOt27WJ0rxY8fnkHXRVT/J6CX+RHW76H2b+B/Rux3e/k7ahbeWrmNmpHlvDvW7szoH19pysUqRQKfpGSIvhmPPzwGtRpxr7hM7jv+xgWb9jCpR0bMOG6ztTT+9RKAFHwS3Db+oO7l79/A3S/g08b3sO4GRspcxXw9PWdGZHSTGOaEnAU/BKcSorg2yfh+1ehTjMOp87ij+lxfLggm/Ob1+GFEV1pGR/tdJUiVULBL8Fn6xL48NfuXX7K7Sxp+1senJXD7sKdPDgoibEDNKYpgU3BL8Gj9Ch88zf3Lj+2GSU3fcjzuQ2Z9HYmLetFM+OeCzm/eZzTVYpUOQW/BIdtS927/H25kHIbOcm/5/6ZuazduZFRPZvzxys6EBWuPw4SHPSbLoGt9Oh/e/m1m+C6eTZv7WrBhDfSiYkI41+3pDCoYwOnqxSpVgp+CVzblnl2+TlwwRh29/ojj8zZwIKcLAa2r8+E67uQEKMxTQk+Cn4JPKXFnl3+K1C7CYz+kE+OtOfxf66gpMzFU9d2ZmQPjWlK8FLwS2DJS3Pv8vdmwwW3cqjf//DE3K3MXLGC5GZ1eGFEMq0TajldpYijFPwSGEqL4bunYPHLENMYRs9iaUhXHnwtnZ0Hj3L/wETuu6QtNTSmKaLglwCQt9yzy18P3X5FycC/8sKCXbw+73ua143i/Xt6c0ELjWmK/EjBL/6rtBi++zssfsm9y795Jrm1e/DA/6WzZschUrs3409XdiQ6Qr/mIhXpT4T4p+3L3dfYyV8H3W7BXjqed1YW8NRbC4mOCGPS6AsY3Kmh01WK+CQFv/iXsmPuXf6iiRDTCG7+gD31+/C7KauYl53PgHYJPD28C/VjajpdqYjPUvCL/6i4yz9/NFz2JJ/nFvHYi/M5WlrO+GvO4+aezTWmKXIaCn7xfWXHYN7TsPBFqNUAbppBYbOL+ctHWcxYnkfnJrG8mNqVNhrTFPGKgl982/YVnl3+Wuh6M1z2JGm7XTz40gK2HzjKvQPa8sCgRI1pipwBBb/4phN3+aPep7TNICZ+lcM/v8ulSVwk0+++kJSWdZ2uVMTvKPjF9+xY6d7l78k6vsvfcDiMB19bzKq8g9xwQVP+fFVHYmrWcLpSEb+k4BffUVYC85+BBf+AWvVh1PvYxEv5z5KtPPlJFpE1Qnn95m4MOa+R05WK+DUFv/iGHemeXf4aSB4FQ55iT1kkj761jG/X59MvKYHnhnehfm2NaYqcKwW/OKusBOY/Cwueh+gEGDUdki5j7ppdPDYzjSPHyvjL1Z245cIWGtMUqSRejUIYY4YYY9YbY3KNMeNOcn+cMWaWMWaVMWapMeY8b9dKENuRDpMudrd3uoyAsT9wpMVAHp2xirvfXU6j2Jp8cn8fftW7pUJfpBKddsdvjAkFXgUuBfKAZcaYOdbarAqHPQ6kW2uvNca09xw/0Mu1EmzKSmDBc+5dflQ8jJwK7YayfMsBHpq+gK37i/jNxW347aAkwsM0pilS2bxp9fQAcq21GwGMMVOBYUDF8O4I/B3AWrvOGNPSGNMAaO3FWgkmO1e5e/m7M6FLKgydQGl4LC9/mc0r3+TQKDaSaXddSI9WGtMUqSreBH8TYFuF23lAzxOOyQCuAxYaY3oALYCmXq4FwBhzF3AXQPPmzb2pXfxJWYl7h7/gOYiqd3yXvzH/MA9OW0xG3kGu69aEJ67uRG2NaYpUKW+C/2TNVXvC7QnARGNMOpAJrATKvFzr/qK1k4BJACkpKSc9RvzUrkyY9WvPLv9GGDIBGxnH5CVb+NvHawkPC+GfN3Xj8s4a0xSpDt4Efx7QrMLtpsCOigdYaw8BYwCM+1m4TZ6PqNOtlQBWXure5c9/1r3LT50C7S8nv/AY495O4+t1e+ibGM+zw5NpGKsxTZHq4k3wLwMSjTGtgO1AKjCq4gHGmDpAkbW2BLgDmG+tPWSMOe1aCVC7Mt3virUrEzqPgKFPQ1RdvsrazaMfrKLwWBl/vrIjt/ZuSUiIJnZEqtNpg99aW2aMuReYC4QCb1pr1xhj7vHc/zrQAXjHGFOO+4nb20+1tmoeiviE8lL3K2/nPwORdeHG96DDlRSVlDF+ZiZTlm6lQ6PaTEntSlKDGKerFQlKxlrfa6enpKTYtLQ0p8uQM7VrtWeXvwo63wBDn4GouqRvK+DBaels3neEu/q15qFLk4gIC3W6WpGAYoxZbq1N8eZYvXJXzl15KSx8AeY9A5F1ju/yy8pdvPpVDi99k0PD2jWZcmcverWu53S1IkFPwS/nZvca9y5/ZwacNxwufxai6rJ57xEenJ7Oyq0FXNO1MX8Zdh6xkRrTFPEFCn45O+VlsOgF+O5p9y5/xLvQ8WqstUxbupW/fpxFWIjh5ZHnc1VyY6erFZEKFPxy5nZneXb56dDpOrj8OYiux77Dx3j0g0y+Wrub3m3q8fyIZBrFRjpdrYicQMEvZ+b7V+GrJyCiNox4BzoOA+DbdXv43YxVHDpayh+v6MBtF7XSmKaIj1Lwi/fS3oS5j0P7K+GqiRAdz9GScp78NIv//LCV9g1j+M8dPWjfsLbTlYrIKSj4xTs5X8Enj0DiYLjhbQgNY1VeAb+dms7GvUe4s28rHh7cjpo1NKYp4usU/HJ6u1bD+7dCg44w/E3KCOG1r3OY+HUOCTERTL6jJ73bxjtdpYh4ScEvp3ZoJ0weARExMGo6Ww6H8OC071mxtYCrkhvzt2HnERulMU0Rf6Lgl1927LA79IsPYsd8yvvry/nLRwsICTFMTO3KsK5NnK5QRM6Cgl9OzlUOH9wOu1dTeN1/eOSrUuauWUWv1nV5fkRXmtTRmKaIv1Lwy8l9/hhkf052yhPcNCeCgqI9PH55e+7o01pjmiJ+TsEvP/fDa7D0DRYnpDJqYRJJDWrw9pgedGysMU2RQKDgl59a9wn288dYGNaLW7ZdyW0XteL3QzSmKRJIFPxyXHneClzTbyPL1Zo/mPt59/ae9EnUmKZIoFHwCwA7Nq8n8p3rOFJeiyltnmHODf2oExXudFkiUgUU/EHOWsvsJevo+NkN1OIYay+Zwt/79cf91skiEohCnC5AnHPgSAn3vbeUup/cSWuzg2PXvcWl/S9W6IsEOAV/kJqfnc9lL8yjb/bf6ReaScjVE0lIvszpskSkGqjVE2SKS8uZ8Nk63lq8mT/EzuXGkG+h7yOEdBvtdGkiUk0U/EFk9faDPDgtnZw9h3m2wwZu2PS2++0SL/mj06WJSDVS8AeBcpdl0vyN/OPL9cRFhTPrqjDO/+Zv0KwXDHsV1NMXCSoK/gCXd6CIh6ZnsHTTfoZ0asjTA2oRO3koxDaB1MlQo6bTJYpINVPwB7CCohJufOMHDh4t5dnhXRjeMRrzf4PBuuCmGRBdz+kSRcQBCv4AZa3lkfcz2FNYzIx7epPcKBLevRYKtsAtc6BeG6dLFBGHaJwzQP3fwk18tXYPj1/egeSmsTDnPtiyCK55DVpc6HR5IuIgBX8AWrH1ABM+W8dlnRpwa++W8N0EWDXNPb3TebjT5YmIwxT8AaagqIT7Jq+kYWxNnhmejMmYCvMmQNeboO8jTpcnIj5APf4A4u7rrzre14/d9YO7xdOqH1z5osY2RQTQjj+guPv6u3lsaAeSa+6BaTdB3dYw4l0I05U2RcRNO/4AsdLT1x/csQFjukbDvwZBaDjc9D5E1nG6PBHxIQr+AHCwqJR7PX39Z4clYaZeB4d3w62fQlwLp8sTER+j4Pdz1loemeGe13//7l7Ezr0X8tJgxDvQ9AKnyxMRH6Qev597c9FmvszazbihHei6/kXImg2Dx0PHq50uTUR8lILfj6VvK2DCZ2u5tGMDbqv5HSyaCCm3w4X3Ol2aiPgwBb+fcvf1V1A/piYvXLAX88nD0PZSGPqMxjZF5JTU4/dD1lp+NyODXQeL+WhEHWrNvgXqd4Qb/g2h+l8qIqemlPBD/160mS+ydvPUoHp0+OY2iKgFo6ZBRIzTpYmIH1Dw+5mMbQX8/bO1XNEuhpEbfg9HC+C2z9zX1xcR8YJXPX5jzBBjzHpjTK4xZtxJ7o81xnxkjMkwxqwxxoypcN9mY0ymMSbdGJNWmcUHm4NHSxk7eQUNa9XgxbBXMLsy4Ya3oFGy06WJiB857Y7fGBMKvApcCuQBy4wxc6y1WRUOGwtkWWuvMsYkAOuNMe9Za0s89w+w1u6t7OKDibWW33v6+ou7zqVG1ly4/DlIGux0aSLiZ7zZ8fcAcq21Gz1BPhUYdsIxFogxxhigFrAfKKvUSoPc24s3M3fNbt7ptIL6WW9Br7HQ406nyxIRP+RN8DcBtlW4nef5WkWvAB2AHUAm8IC11uW5zwJfGGOWG2Pu+qUfYoy5yxiTZoxJy8/P9/oBBINVeQU8+elaHm6ey4U5z0H7K90v0hIROQveBP/JhsLtCbcvA9KBxkBX4BVjTG3PfRdZa7sBQ4Gxxph+J/sh1tpJ1toUa21KQkKCN7UHhR/7+n2i8rj3wARM4/Phuv+FkFCnSxMRP+VN8OcBzSrcbop7Z1/RGGCmdcsFNgHtAay1Ozz/3QPMwt06Ei9Yaxn3wSooyGNS2DOYqHgYORXCo5wuTUT8mDfBvwxINMa0MsaEA6nAnBOO2QoMBDDGNADaARuNMdHGmBjP16OBwcDqyio+0L3z/RYWrt7I7LiJ1HCVwE3TIaaB02WJiJ877VSPtbbMGHMvMBcIBd601q4xxtzjuf91YDzwljEmE3dr6FFr7V5jTGtglvs5X8KAydbaz6vosQSUzLyDPP1JJtPrvEbc0c1w0wyo38HpskQkAHj1Ai5r7afApyd87fUKn+/AvZs/cd1GQEPmZ+hQcSlj31vOUxFvc17xcrj6FWgzwOmyRCRA6CJtPubHvv4VhdO5xvUl9H0Yuo12uiwRCSAKfh/z7g9bYM2HPBo2BTpdBwP+6HRJIhJgdK0eH7J6+0E+/mQ274W/hm3WC3PNaxCiv5tFpHIpVXzEoeJS/vafT3kj7DlC6jTBpE6GGjWdLktEApCC3wdYaxk/fRFPFv2VmIgQQm/+AKLrOV2WiAQotXp8wOTFOVyfO46WofmEjpoD8W2dLklEApiC32Gr8wqInvsQvULW4rrmf6FFb6dLEpEAp1aPgwqLS1n69qNcE7KAoovGEZI8wumSRCQIKPgdYq1l5lvPc1vpVPLbDidq0M/e30ZEpEoo+B3y5WczGbnzGfLqpJCQ+hqYk10EVUSk8in4HZC9ZgU9ltzH3hqNaXznDAgLd7okEQkiCv5qdnj/TqJnpOIyYUTeOpOQ6DinSxKRIKPgr0a2pIg9k66nnms/Oy//N3FNk5wuSUSCkIK/urhcbH3zV7Q8msW3nZ6kU4+BTlckIkFKwV9N8mc/TotdXzAt7k4uG/6Lbz0sIlLlFPzVoHjJmyRkvMYHIYO59PbxhIRogkdEnKPgr2I292tqfPYw37mSaTrqFeJjdOE1EXGWgr8q7V5D2ZTRZLuasr7PRHq21fvliojzFPxVpXAXpe8MZ39ZOK83eYo7B3V1uiIREUAXaasaJUcof28EZUf283CN8bx402D19UXEZ2jHX9lc5dgZt2F2ZXJf6b3cO+p64mtFOF2ViMhxCv7KNvcPmOzPeaJ0NMmXpNKrtd5QRUR8i4K/Mi15A5a8xluuoWxqfRO/GaA3VBER36Mef2VZ/xn283EsCu3J62Fj+PjGroSqry8iPkg7/sqwYyV2xm1si0jk7qK7eWFkivr6IuKzFPzn6mAeTE6lKCyW6wse4O5BXbiwjfr6IuK7FPznovgQvDeC8pIj3Fj4EO3atmWs+voi4uPU4z9b5aXw/q+we9czLuJP7I5szb/V1xcRP6Ad/9mwFj59BDZ8w7QGD/FBQSITU7uSEKO+voj4PgX/2Vg0EZa/xdo2dzBuU1ceGJhE7zbxTlclIuIVBf+ZWvMhfPU/HGpzNddlX0LvNvW49xL19UXEf6jHfya2LYNZd1PepAcjdo8mOiKEF1PV1xcR/6Idv7f2b4IpqdiYRjwR9QfW7yvlpdSu1Nf19UXEzyj4vXH0ALx3A7jK+Dz5Jd7NPML9lyTSu636+iLifxT8p1NWAtNGQ8EWtl32Lx78+ggXtq7H/QMTna5MROSsqMd/KtbCR/fD5gUcu/p1xnwbTq0Iw8SR6uuLiP/Sjv9U5j0DGVPg4sd5PLcjG/IPM1F9fRHxcwr+X5IxDb57CpJH8n70SD5Ykcd9lyRykfr6IuLnvAp+Y8wQY8x6Y0yuMWbcSe6PNcZ8ZIzJMMasMcaM8XatT9q8CGaPhZZ9ye75JH+as4ZerevygPr6IhIAThv8xphQ4FVgKNARGGmM6XjCYWOBLGttMnAx8LwxJtzLtb5lbw5MHQV1W1F07VuMnbqaWhFhvJR6vvr6IhIQvNnx9wByrbUbrbUlwFRg2AnHWCDGGGOAWsB+oMzLtb7jyF54bziEhMGo6fz5i+3k5h/mhRu7Ur+2+voiEhi8Cf4mwLYKt/M8X6voFaADsAPIBB6w1rq8XAuAMeYuY0yaMSYtPz/fy/IrUWkxTBkJhbtg5FRmbKrBjOV53DegLX0TE6q/HhGRKuJN8J+sv2FPuH0ZkA40BroCrxhjanu51v1FaydZa1OstSkJCQ4E7bdPQt5SuPYNcsLb86cPV9OzVV0eGJRU/bWIiFQhb+b484BmFW43xb2zr2gMMMFaa4FcY8wmoL2Xa31D34eg8fkcTbyKsa8uJCo8lJdGqq8vIoHHmx3/MiDRGNPKGBMOpAJzTjhmKzAQwBjTAGgHbPRyrW+IjIPzruN/5qwmZ4+7r99AfX0RCUCn3fFba8uMMfcCc4FQ4E1r7RpjzD2e+18HxgNvGWMycbd3HrXW7gU42dqqeSjn7oPleUxPy+O+S9rSL0l9fREJTMbdnfEtKSkpNi0trVp/Zu6eQq56eRGdm8Yy+Y6ehIXqtW0i4j+MMcuttSneHKt0A46WlDP2vZVEhYfy8sjzFfoiEtB0kTbgiTlryN5TyNtjeqivLyIBL+i3trNW5jEtbRu/ubiN+voiEhSCOvhz9xzmD7NW06NlXR7UvL6IBImgDX53X38FNWu45/XV1xeRYBG0Pf6/fLSG9bsLeWtMdxrGqq8vIsEjKLe5H67cztRl7r7+xe3qO12OiEi1Crrg35B/mMdnZdK9ZRwPXaq+vogEn6AK/uJSd18/IixEfX0RCVpB1eP/y0drWLerkH+P6U6j2EinyxERcUTQbHlnp29nytJt/PriNgxQX19EglhQBP+G/MM8PjOTlBZxPKy+vogEuYAP/h/7+uHq64uIAEHQ4//rx1nuvv6t3WlcR319EZGA3v7OTt/O5CVbubt/awa0V19fRAQCOPg3evr6F7SI45HB7ZwuR0TEZwRk8BeXljN28kpqhIXw8sjzqaG+vojIcQHZ4x//cRZrdx7izVtT1NcXETlBwG2FP8rYwXtLtnJ3v9Zc0r6B0+WIiPicgAr+TXuP8NjMTLo1r8Mjl6mvLyJyMgET/D/O64eFGl4e1U19fRGRXxAwPX6XtbRvFMPDg5Noor6+iMgvCpjgjwoP4x8jujpdhoiIz1M/REQkyCj4RUSCjIJfRCTIKPhFRIKMgl9EJMgo+EVEgoyCX0QkyCj4RUSCjLHWOl3Dzxhj8oEtZ7k8HthbieVUNX+rF1RzdfG3mv2tXgismltYaxO8+QY+GfznwhiTZq1NcboOb/lbvaCaq4u/1exv9ULw1qxWj4hIkFHwi4gEmUAM/klOF3CG/K1eUM3Vxd9q9rd6IUhrDrgev4iInFog7vhFROQUFPwiIkHGb4LfGDPEGLPeGJNrjBl3kvuNMeYlz/2rjDHdvF3rozVvNsZkGmPSjTFpPlRze2PM98aYY8aYR85krQ/W66vn+CbP78MqY8xiY0yyt2t9tOZqP89e1DvMU2u6MSbNGNPH27U+WvOZnWNrrc9/AKHABqA1EA5kAB1POOZy4DPAAL2AJd6u9bWaPfdtBuJ98DzXB7oDTwKPnMlaX6rXx89xbyDO8/lQP/ldPmnNTpxnL+utxX+f4+wCrPODc3zSms/mHPvLjr8HkGut3WitLQGmAsNOOGYY8I51+wGoY4xp5OVaX6vZKaet2Vq7x1q7DCg907U+Vq9TvKl5sbX2gOfmD0BTb9f6YM1O8Kbew9aTmEA0YL1d64M1nzF/Cf4mwLYKt/M8X/PmGG/WVoVzqRnc/1O/MMYsN8bcVWVVel9PVa49W+f6M/3hHN+O+1+FZ7O2spxLzVD959mreo0x1xpj1gGfALedydoqcC41wxmeY395s3Vzkq+d+LfdLx3jzdqqcC41A1xkrd1hjKkPfGmMWWetnV+pFf7cuZwrJ87zuf5Mnz7HxpgBuEP0x16uL/8uuw/8ec1Q/efZq3qttbOAWcaYfsB4YJC3a6vAudQMZ3iO/WXHnwc0q3C7KbDDy2O8WVsVzqVmrLU//ncPMAv3PwWr2rmcKyfO8zn9TF8+x8aYLsC/gGHW2n1nsrYKnEvNTpznMzpPnoBsY4yJP9O1lehcaj7zc1zVT1pUxgfuf5lsBFrx3yc+Op1wzBX89InSpd6u9cGao4GYCp8vBob4Qs0Vjn2Cnz65W+3n+Rzr9dlzDDQHcoHeZ/t4fajmaj/PXtbblv8+UdoN2O75c+jL5/iXaj7jc1ylD6aST8zlQDbuZ77/4PnaPcA9ns8N8Krn/kwg5VRrfblm3M/sZ3g+1vhYzQ1x704OAQWez2s7dZ7Ptl4fP8f/Ag4A6Z6PND/4XT5pzU6dZy/qfdRTTzrwPdDHD87xSWs+m3OsSzaIiAQZf+nxi4hIJVHwi4gEGQW/iEiQUfCLiAQZBb+ISJBR8IuIBBkFv4hIkPl/L/MrAMokbBgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hs_eu = [prod(alphas_us,m) for m in ms_eu]\n",
    "ms_eu = ms_eu[:]\n",
    "hs_eu = hs_eu[:]\n",
    "plt.figure()\n",
    "plt.plot(ms_us,hs_us)\n",
    "plt.plot(ms_eu,hs_eu)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1b2eb0f",
   "metadata": {},
   "source": [
    "## Demand : Price Elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d9341d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us*1.01,mu_us,s_us) for i in yy_us]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b28236ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9.2007609786807, -0.7265809284720925)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share0 = p_us*np.mean(ms0)/np.mean(yy_us)\n",
    "share1 = p_us*1.1*np.mean(ms1)/np.mean(yy_us)\n",
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (p_us*1.01-p_us)/p_us\n",
    "\n",
    "((share1 - share0)/share0)/dxox, ((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "786586f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us*2,mu_us,s_us) for i in yy_us]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us*2*1.01,mu_us,s_us) for i in yy_us]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "e6842fc2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9.432891805429465, -0.5155529041550285)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share0 = p_us*2*np.mean(ms0)/np.mean(yy_us)\n",
    "share1 = p_us*2*1.1*np.mean(ms1)/np.mean(yy_us)\n",
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (p_us*2*1.01-p_us*2)/(p_us*2)\n",
    "\n",
    "((share1 - share0)/share0)/dxox, ((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "305d34c7",
   "metadata": {},
   "source": [
    "## Demand: Rand Experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "7ef9823e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def focr(m,sigma_c,phi,alphas,y,p,mu,mom,mu0):\n",
    "\treturn mu*p*( (1-(1-mu0)*mom)*y - mu*p*m )**(-sigma_c) - alphas[0]*phi*np.exp(-alphas[1]-alphas[0]*m)\n",
    "\n",
    "def ur(m,sigma_c,phi,alphas,y,p,mu,mom,mu0):\n",
    "\treturn ( ((1-(1-mu0)*mom)*y - mu*p*m)**(1-sigma_c) )/(1-sigma_c) + phi*(1-np.exp(-alphas[1] - alphas[0]*m))\n",
    "\n",
    "def optmr(sigma_c,phi,alphas,y,p,mu,mom,mu0):\n",
    "\ttry:\n",
    "\t\tm = bisect(focr,0.0,y*0.99,args=(sigma_c,phi,alphas,y,p,mu,mom,mu0))\n",
    "\texcept:\n",
    "\t\tup  = ur(0.99*y,sigma_c,phi,alphas,y,p,mu,mom,mu0)\n",
    "\t\tlow = ur(0.0,sigma_c,phi,alphas,y,p,mu,mom,mu0)\n",
    "\t\tif up > low:\n",
    "\t\t\tm = 0.99*y\n",
    "\t\telse :\n",
    "\t\t\tm = 0.0\n",
    "\treturn m\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "eda6d540",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optmr(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us,mu_us) for i in yy_us]\n",
    "ms1 = [optmr(sigma_c,phi,alphas_us,i,p_us,mu_us*1.01,s_us,mu_us) for i in yy_us]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "8588c849",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.7265809284720925, -0.7265809284720925)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share0 = p_us*np.mean(ms0)/np.mean(yy_us)\n",
    "share1 = p_us*np.mean(ms1)/np.mean(yy_us)\n",
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (mu_us*1.01-mu_us)/mu_us\n",
    "\n",
    "((share1 - share0)/share0)/dxox, ((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ee1b5c3",
   "metadata": {},
   "source": [
    "## Demand : Income Elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "4198d7fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us*1.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "892c5f2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2541301907224377, 1.2566714926296656)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share0 = p_us*np.mean(ms0)/np.mean(yy_us)\n",
    "share1 = p_us*np.mean(ms1)/np.mean(1.01*yy_us)\n",
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (np.mean(1.01*yy_us)-np.mean(yy_us))/np.mean(yy_us)\n",
    "\n",
    "((share1 - share0)/share0)/dxox, ((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "9f47ac85",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us*2]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us*2*1.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "c4dc5d9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.14451965456975438, 0.8540351488845399)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share0 = p_us*np.mean(ms0)/np.mean(yy_us*2)\n",
    "share1 = p_us*np.mean(ms1)/np.mean(1.01*yy_us*2)\n",
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (np.mean(1.01*yy_us*2)-np.mean(yy_us*2))/np.mean(yy_us*2)\n",
    "\n",
    "((share1 - share0)/share0)/dxox, ((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "183a75b8",
   "metadata": {},
   "source": [
    "## Demand : Price Elasticity in Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "f10b4443",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu]\n",
    "ms1 = [optm(sigma_c,phi,alphas_eu,i,p_eu*1.01,mu_eu,s_eu) for i in yy_eu]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2d2c29ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8.939059782513258, -0.9644911068061123)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share0 = p_eu*np.mean(ms0)/np.mean(yy_eu)\n",
    "share1 = p_eu*1.1*np.mean(ms1)/np.mean(yy_eu)\n",
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (p_eu*1.01-p_eu)/p_eu\n",
    "\n",
    "((share1 - share0)/share0)/dxox, ((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2469bccd",
   "metadata": {},
   "source": [
    "## Demand : Income Elasticity in Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "02f0e64d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu]\n",
    "ms1 = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu*1.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "744a8c18",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.7700534600029438, 1.7777539946029677)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "share0 = p_eu*np.mean(ms0)/np.mean(yy_eu)\n",
    "share1 = p_eu*np.mean(ms1)/np.mean(1.01*yy_eu)\n",
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (np.mean(1.01*yy_eu)-np.mean(yy_eu))/np.mean(yy_eu)\n",
    "\n",
    "((share1 - share0)/share0)/dxox, ((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3762e4c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "cf2a50979671a58939829e6829efb726aa5da11149213b77bd50351f899d04fb"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
